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ABSTRACT 

The  partial  differential  equations  for  supersonic  conical  flow 
are  nonlinear  and  hyperbolic.     By  inverting  the  problem  and  assuming 
the  shape  of  the  shock,   there  results  a  Cauchy  problem  with  initial 
data  given  by  the  Rankine-Hugoniot  conditions.     This  Cauchy  problem 
is  then  solved  numerically  by  means  of  a  program  written  for  the 
IBM  7090  Computer.     The  numerical  method  was  developed  by  R.    D. 
Richtmyer  and  G.    E.    Lewis.     Customarily,    finite  difference  schemes 
are  tried  to  solve  these  types  of  problems.     However,    Richtmyer  and 
Lewis  utilized  a  power  series  technique  suggested  by  the  Cauchy- 
Kowalewski  theory.     The  numerical  procedures  involve  algorithms 
which  direct  the  machine  to  compute  the  power  series,    techniques  of 
significance  arithmetic,    and  methods  of  analytic  continuation.     The 
feature  of  the  methods  involved  is  the  high  degree  of  accuracy  in  the 
computation. 
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NYO-10,  432 

A  POWER  SERIES  SOLUTION  FOR  COMPRESSIBLE  FLOW 
PAST  A  CONICAL  SHOCK  WAVE 

by  ■ 

Jay  A.    Leavitt 

1.     Introduction  . 

The  usefulness  of  power  series  for  solving  differential  equa- 
tions has  been  recognized  for  more  than  a  century.     The  Cauchy- 
Kowalewski  Theorem  proves  the  existence  of  a  unique  analytic  solution 
for  noncharacteristic  initial  value  problems  and  at  the  same  time  pro- 
vides a  method  for  finding  a  power  series  solution.     However,    until 
recently,   this  method  has  been  useless  as  a  computational  technique 
except  for  the  simplest  of  problems.     In  order  to  apply  the  Cauchy- 
Kowalewski  Theorem  to  a  supersonic  flow  involving  a  shock,    it  is  first 
necessary  to  describe  the  flow  by  the  "inverse  method".     In  the     in- 
verse method"  one  assumes  the  shape  of  the  shock  and  the  conditions 
ahead  of  it,    and  tries  to  find  the  flow  behind  the  shock.     By  appropriate 
choice  of  coordinates  the  shock  acts  as  the  initial  curve,     ?  =  0;      the 
Rankine-Hugoniot  conditions  provide  the  initial  conditions  for  the  sys- 
tem of  equations  describing  the  flow.     The  Cauchy-Kowalewski  method 
is  a  recursive  technique  for  finding  the  higher  order  derivative  of  the 
unknowns  in  terms  of  the  lower  order  ones.     Several  different  attempts 
at  solving  fluid  dynamics  problems  by  firjding  truncated  power  series 
failed  because  of  singularities  of  the  solutions.     These  methods  failed 
to  succeed  whenever  the  singularities  were  closer  to  the  point  of  ex- 
pansion than  the  body.     These  difficulties  were  discussed  by  Van  Dyke 
[1].       In  his  paper,   he  states  that  for  the  case  of  axisymmetric  flow  past 
a  blunt  body,    there  exists  a  limiting  line,    a  line  which  is  an  envelope  of 


characteristics.     The  first  extensive  success  using  power  series  meth- 
ods was  initiated  by  R.    D.    Richtmyer  on  the  Univac  and  was  later  com- 
pleted by  G.    Lewis  on  the  I.  B.  M.    704. [2,    3].       In  addition  to  providing 
the  recursive  relations  by  which  the  machine  could  calculate  the  series 
coefficients  by  subroutines  Richtmyer's  main  idea  was  to  carry  an 
accuracy  index  with  each  coefficient.     Those  terms  where  significance 
was  lost  were  discarded  so  that  they  would  not  interfere  with  the  calcu- 
lation.    This  was  done  by  inspecting  the  accuracy  index  as  each  term 
was  calculated.     G.    Lewis  was  able  to  complete  this  problem  by  devis- 
ing two  methods  of  continuation  for  truncated  power  series.     The  diffi- 
culty lay  in  the  fact  that  a  truncated  power  series  is  only  a  polynomial, 
and  to  expand  it  around  another  point  is  just  a  rearrangement  of  terms, 
which  provides  nothing  new.     The  proofs  of  the  convergence  of  his 
methods  of  analytic  continuation  are  found  in  his  thesis  [3]. 

After  the  bow  wave  problem  was  solved,   G.   Whitham  sug- 
gested that  the  method  be  applied  to  the  conical  flow  problem.     It  is  to 
this  end  that  this  paper  is  directed.     Solutions  to  a  linearized  version 
of  this  problem  were  presented  by  A.    Ferri  [4].      Another  technique  was 
used  by  B.    Briggs  [5,    6],    who  described  the  flow  in  terms  of  two  stream 
functions.     He  used  a  difference  method  to  integrate  the  resulting  equa- 
tions. 

The  author  wishes  to  thank  R.    D.   Richtmyer  whose  encourage- 
ment and  guidance  made  this  paper  possible. 

2,     Methods  of  Continuation 

In  this  problem  both  of  Lewis'  methods  were  tried.     The  second 
was  used  in  the  final  computation.     The  first  method  was  discarded  be- 
cause it  was  found  that  in  repeated  continuations,   the  coefficients  lost 
most  of  their  significance.     In  a  test  where  a  series  was  expanded  about 
a  point    x  ,     and  then  continued  as  an  expansion  about    x. ,     it  was  found 


that  when  this  second  series  was  re -expanded  around    x   ,     about  half  of 

o 

the  terms  were  lost  due  to  the  loss  of  significance.     The  second  method 
proved  quite  satisfactory*   especially  since  it  required  no  new  coding. 
The  program  that  computed  the  power  series  from  the  initial  data  was 
simply  reapplied  with  new  initial  data  obtained  from  the  series  to  be  con- 
tinued.    The  second  method  has  the  additional  advantage  of  making  more 
use  of  the  differential  equations  than  the  first  method. 

a.     First  method  of  continuation 

Let    F(x)    be  a  function  of    x  -  (x^ ,  .  .  . ,  x   ),     analytic  in  each 
variable,    which  is  expressible  as  a  power  series  about  the  origin.     Let 
D      be  the  domain  of  convergence  of  this  power  series  and  let    x      be 
any  fixed  point  in    D  .       F(x)  may,   then,   be  expressed  by  a  power  series 

about    x  .      Let    D      be  the  domain  of  convergence  of  this  new  power 

1 
series  and  let    x      be  any  fixed  point  in    D  .      The  problem  is  to  find  an 

1 
approximation  to    F(x   )    when  only  the  truncated  power  series  about  the 

origin  of    F{x)    is  known. 

Introduce  the  new  coordinate  system 

o  ,^  o 


and  let 


X     —     X 


=  o  1  o 

;       =    X      —    X 


(C.   =  X.  -  x") 

111 


Let    f(?)  =  F(?  +  x°),     then    F{x^)  =  f(?°).       Denote  by    F   (x)    the  known 

truncation  of  the  expansion  of    F(x)    about  the  origin  to  order    k    and  let 

f,  (?)  =  F,  (C+x°).       Denote  by    f     ,  (?)    those  terms  of    f,  (?)    of  order 
k  k  s,  k  k 

less  than  or  equal  to  [sk]  when  s  is  some  constant  between  zero  and 
one  ([sk]  is  the  greatest  integer  in  sk).  This  is  a  further  truncation 
of    F   (x)    written  as  an  expansion  of    ?. 


Theorem:      Under  these  assumptions,    there  exists  an    s, 
0   <   s  <    1,     such  that 

lim    f     ,  (C°)    =    F(x^)      . 

s,  k 
k—  00 

The  theorem  shows  that  an  analytic  continuation  can  be  found 

since    x     e   D,     is  not  necessarily  in    D   . 
1  o 


b.     Second  method  of  continuation 
Let 

(A)  1^     =    G(u)  1^         and  u(0,  y)    =    f(y)     , 

9x  9y 

°° 

where    G(u)    is  analytic  at    u  =  0,     G(y)  =    X       a    u   .       (In  the  problem, 

n=0       ^ 
G    also  depends  explicitly  on    x    and    y,     but  for  simplicity  assume 

G  =  G(u).       The  theorem  is  true  for  the  most  general  analytic  non- 
characteristic  initial  value  problem.  )     f(y)    is  analytic  at    y  =  0, 

00 


Z-  " 


^(y^  ^    /         t)    y      and  only  the  coefficients  of  order  less  than  or  equal 

n  =  0       " 

to    k    of  this  expansion  of    f(y)    are  known.     Let    u   (x,  y)    be  the  poly- 
nomial of  order    k    obtained  from  the  series  expansion  of    u(x,  y)    in 
powers  of    x    and    y. 

Let    (x  ,  y   )    be  a  point  in  the  domain    D      of  absolute  conver-- 
o  -^  o  ^  o 

gence  of  the  power  series  of    u(x,  y)    expanded  about  the  origin,    and  let 
(x  ,  y   )    be  a  point  in  the  domain    D      of  absolute  convergence  of  the 
power  series  of    u(x,  y)    expanded  about  the  point    (x   ,y   ).       The  prob- 
lem is  to  find  an  approximation  to    u(x  ,  y   )    when  only 

k 


yy)  =    Zb„y"    , 


n  =  0 
and  therefore  only    u   (x,  y),     is  known.     As  in  the  first  method 


{x.,y^)   €   D      need  not  be  in    D   .       The  procedure  for  finding    u   (x,  y) 

111  O  K 

is  found  in  Section  7. 


Introduce  the  new  coordinates 

I   =  X   -  x^    ,  ?    =  y   -  y^ 

^o   =  ^    -  ^o    '  ^o   =  ^1    -  ^o 

and  let 

If    =    G(w)  If     ,  w(0,  ?)     =    u,(x.?+y^)      .  (B) 

OZ  05  K       O  O 

Also  let    w     ,(?,?)    be  that  polynomial  which  is  obtained  from  the  series 
expansion  of    w(C,  ?)    in  powers  of    |    and    g    where  all  the  terms  of 
order  greater  than    [sk]    are  discarded,     s    is  some  constant  between 
zero  and  one.     G.    Lewis  found  that    s  =  1     was  apparently  best  for  most 
cases,    but  was  only  able  to  prove  this  theorem  for    s   e    (0,  1).       The 
initial  value  problem  (B)  is  in  the  same  form  as  (A).     The  only  differ- 
ence is  that  the  initial  data  is  given  on  the  line    x  =  x   ,       The  coeffi- 

°  o 

cients  to  order    k    of    w(0,  ?)    are  known  from    u  (x  ,  ?+y   ).       The  same 

K.       O  O 

algorithms  that  are  used  to  find  the  coefficients  of    u   (x,  y)    can  be  used 

K. 

to  find  the  coefficients  for  the  series  of    w(C,  0= 
Finally  consider  the  problem 

|J     =    G(v)  ly     ,  v(0,y)    -•  l(y)      ,  (C) 

00  00 

where    G(v)  -    ^        |a    |  v       and    Uy)  =    /^        |b    |  y   . 
n-0  ^  n=0  ^ 

Theorem:      If    (x  ,  y   )    is  in  the  domain  of  convergence  of  the 

power  series  of    v(x,  y),     the  solution  to  problem  (C),    expanded  about 

the  origin;     and  if    (x  ,  y^  )    is  in  the  domain  of  convergence  of  the  power 

series  of    u(x,  y)    expanded  about  the  point    (x   ,  y   ),     then  there  exists 

an    s   e   (0,  1)    such  that 

lim    w     ,(?,?)=    u(x,,  y,  )      . 

s,  k     o      o  1    "^  1 

k—  00 
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3.     The  Problem 

The  general  equations  of  motion  for  an  inviscid  compressible 
flow  in  spherical  coordinates  are  given  by 

2         2 
9w         u   9w    ,  V        9w    ,     1    9p  (u     +  V    )  „ 

9r  r    90  r  sm  0  9(p  p    9r  r 

2 
9u         u    9u  V        9u     4.    J_   9p     ,     (uw  -  v     cot  0)    _ 

^   9r  r    90  r  sin  0  9(|)  pr    90  r 

(1) 

9v         u    9v  V        9v    ^  1         9p_    ^    (vw  -^  uv  cot  0)      _ 

9r  r    90  r  sin  9  d(j)  pr  sin  Q  d^  r 

•;r—  (pr     w  sin  0)    +    ^tt  ^pru  sin  0)    +    -r-r  (prv)      =    0    , 
9r  90  o0 

where    w  -  dr/dt,     u  =  r  d0/dt,     v  -  r  sm  0  d(/)/dt,     p  -  pressure  and 
p  =  density. 

The  first  three  equations  are  the  equations  for  conservation  of 
momentum;    the  fourth  is  the  continuity  equation.     To  these  must  be 
added  the  energy  equation  for  a  steady  adiabatic  flow;    along  a  stream- 
line   u  9S/90  +  (v/sin  0)(9S/90)  =  0,     where    S    is  the  entropy.     The  equa- 
tion of  state  for  a  polytropic  gas  is 

(2)  p    =    A(S)p^      , 

where    y     is  the  gas  constant.      S    is  removed  from  the  energy  equation 
using  (2).     This  leaves  five  differential  equations  in  five  unknowns. 

Inside  the  shock,   the  flow  is  considered  conical,   that  is,   the 
condition  that  ail  flow  quantities  are  independent  of  the  variable    r. 
This  assumption  implies  the  following  two  equivalent  statements.     Along 
each  ray  through    r  =^  0,     the  flow  quantities  are  equal  at  every  point. 
It  is  sufficient  to  solve  the  equations  on  the  spherical  surface    r  -  con- 
stant,    for  the  flow  is  similar  on  all  spherical  surfaces  (actually  identi- 
cal except  for  a  scaling  factor). 


attached  conical  shock 


Figure  1 
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Let  the  shock  be  given  by 

(3)  g(<S>)    -    9      . 

In  order  to  use  the  Cauchy-Kowalewski  method,   the  equations  must  be 
rewritten  in  terms  of  new  coordinates  such  that  the  shock  becomes  the 
initial  curve    C  =  0    (the  shock  is  actually  the  cone    g((/))  =  0,     r  =  r,     but 
we  are  only  considering  the  intersection  of  this  cone  with  the  surface 
r  =  constant    so  that  the  shock  is  represented  by  a  curve). 

Let 

(4)  ?     -    g{<i>)    -    0      ,  r)    =    (j)       . 

The  shock  is  now  the  line    1=0.      The  equations  inside  the  shock  be- 
come: 

|£    .     i_[vf^_T^g,]|P   +^sin0  L|^-v|^+[ucos0  -wsin0]u 
9?  (  p     *=      9n  (      9r]  9rj 

+  [g'w  +  vcos0]v  ~  2f2w|  \/\^^  ~  —  (g'^  +  sin^0)  |     , 
<  -  V  —    +   (u    +  v    )  sin  0  > 


(5) 


9w 

9f    "  n  ' 

9u  0^2  ^   ^    sin  0    9p 

;  -  V   -T—  -  uw  Sin  0   +  V    cos  0   +    — — 

8u    _     ( drj p        9g 

9?     "  ^ 


9v  .      n 

,    -  V  -r—   -  vw  sm  0  —  uv  COS 
9v  (  dri 


0  _   1  9p    _   g^  9p  ) 
p    9n  p    9g  ) 


9C  CI 

9p     ^     ( 9r)  7p  9g dri     ) 

9?     "  n 

where    57  =  vg'  -  u  sin  0,     g'  =  dg/dri,     sin  0  =  sin  (g(n)  -?), 

cos  0  =  cos  igir))  -i),     p  =  p(C,  n),     P  =  P(?.  n),     u  =  u(C,  ri),     v  =  v(g,  rj), 

w  =  w(g,  TTj).      The  characteristic  equations  for  this  system  are 


-  11 


(6) 


+  n(<n>^-   ^(g'^+  sin^0))dr7^  I 
P  ) 


where 


§-     =     ?— ^(vg'  -  u  sin  G)       and       ^     =     ^     .  (7) 

dt  r  sin  0      "^  dt  r  sin  0 

The  form  (5)  of  the  flow  equations  is  valid  so  long  as  neither 

(a)  0=0    nor 

(b)  r2  =  0    nor 

(c)  n^  -  ^(g'^+  sin^0)  =  0. 

In  case  (a),     0  =  0    corresponds  to  the  origin  of  the  polar  co- 
ordinates   (0,  (/))    on  the  spherical  surface  (this  point  will  be  referred  to 
as  the  polar  origin;      0    and  there    C    correspond  to  the  radial  component, 
</)    and    r]    refer  to  the  angular  one).     The  polar  origin  lies  inside  the 
body.     If  either  case  (b)  or  (c)  should  occur,   the  curve    ?  =  constant 
would  be  characteristic,    which  can  be  seen  by  setting    d?  =  0    in  (6). 

This  would  violate  the  conditions  of  the  Cauchy-Kowalewski  Theorem. 

/2  2 

sin  g  +  g'      is  the  normal  com- 
ponent,    N,     of  velocity  (see  below).      N  =  0    only  in  the  case  of  a  con- 
tact discontinuity;    by  the  conservation  of  mass  it  cannot  vanish  other- 
wise.    If  (c)  were  the  case,   then  ■ 


I 


•    2     ,      ,2/  ^2 

sin   g  +  g'  /  _      N 

TP  7P 

P  P 


t 

This  is  a  contradiction  since  the  normal  component  of  velocity  must  be 
subsonic  behind  a  stationary  shock.     The  stagnation  point  on  the  body 
is  an  exceptional  point  for  the  solution,    for  there    N  =  0    and    dC  =  0, 


12 


so  that  it  cannot  be  used  as  a  point  of  continuation  by  the  second  method. 
In  terms  of    u    and    v    the  flow  was  found  to  be  subsonic  between  the 
shock  and  the  body,    so  there  are  probably  no  exceptional  points  interior 
to  the  flow  where  a  continuation  would  not  be  possible. 

Figure  2a  is  a  cross -section  picture  looking  at  the  plane 
^  =  constant.       Figure  2b  is  a  drawing  looking  down  on  the  spherical  sur- 
face   r  =  constant.       From  Figure  2c,    sin  a    and    cos  a    are  readily 
found.     Since    0  =  g(0),     d0  =  g'  d(|) , 

tQ\                      ■                              g'                                                           sin  (g) 
(8)  sin  a    =    s ,  cos  a    -    ^ 


/       2  2  /       2 

\  sin  (g)  +  g'  ysin  (g) 

so  that  the  components  of  velocity  at  the  shock  are 

u  sin  (g)  -  vg' 


N    =     -  [(-u)  cos  a  +  V  sin  a] 


(9) 


,/    .    2,    ,  72 

y  sm   (g)  +  g' 

/       \     •  ,  ug'  +-  V  sin  (g) 

(-u)  sm  a    i-    V  cos  a    =     —  - 


sin   (g)  +  g' 

where    N    is  the  normal  component  of  velocity  and    T    the  tangential. 
The  conditions  ahead  of  the  shock  are 

p^  =    1,       Pq  =    1,       7   =    1.405,       M    =    5.4223,       Q   -    M /7, 

(10)         u    .=    -  Q  sin  0,       V     ^    0,       w     =   Q  cos  0, 
o  o  o 


27       ^o  2  2  2 

B  - —  +    u  +    V  +    w  ,       M    is  the  mach  number. 

T  -  1    P  o  o  o 

o 

Based  on  these,    we  can  find  the  initial  values  of  the  flow  quantities  be- 
hind the  shock  by  the  Rankine-Hugoniot  jump  conditions: 
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Figure  2 
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O 

^  (X  +  1)(1  +  {-^ — )    )  +  u 

sm  g  o 


Pi 


Xp^    -  1 


(11)  '  ^-'1 


A    -    1  O 

U,      =     U        - 


+      ^^^ 


1  °  ^  (IM-S^)')  ^% 

sm  g 

V,     =    (u    -  u.)(-^)     , 
1  o         1     sm  g 

w.     =    w        , 
1  o 

where    X  =  (7  + 1)/ (7  -  1 ),     g'  =dg/d(/).      If    9    and    (/>    are  considered  as 
radial  and  angular  components  of  polar  coordinates,   then    g    was  chosen 
so  that    9  =  g{<j))    would  be  the  polar  representation  of  an  ellipse.     The 
plane  cross-section  of  the  shock  is  given  by 

R    =    r  cos  (a)  tan  (gW)     , 

where    R    is  the  radial  distance  measured  on  the  plane  surface;    the 
plane  is  normal  to  the  axis  of  the  cone  and  passes  through  the  point 
(r=r,  9=g(0)=a,  (^=0).      In  principle,    any  analytic  noncharacteristic  curve 
would  do.     Two  problems  were  solved.     In  one  case  the  ratio  of    b/a, 
the  minor  to  the  major  axis,    was    1    and  in  the  other  case    b/a  =  .7. 
The  curve    g    is  given  by 

(12)  9    =    gW    = 


The  conditions  ahead  were  normalized  to  make  the  Bernoulli  constant 
equal  to  unity  by 


(13) 


Pi 
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P 
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u    =     

V, 
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1 

V     = 

w    ■■ 
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(note    p  =  p(0,  ri),     p  =  p(0,  rj),     etc.  ).     The  particular  choice  of    M,    7 
and  the  normalization  was  taken  in  order  to  compare  the  circular  case 
(b/a  =  1)    with  Kopal  [7].     The  results  of  this  calculation  are  shown  in 
Table  1  in  Section  8. 

All  the  functions  are  considered  as  power  series  in  two  varia- 
bles,    I    and    rj.       The  equations  (5)  establish  the  recursive  relations 
between  the  coefficients,     and  the  conditions  (11)  and  (13),    found  from 
the  Rankine-Hugoniot  jump  conditions,   provide  the  starting  values. 

4.     The  Algorithms 

Subroutines  automatically  carry  out  the  multiplication,    division, 
differentiation,    and  integration.     These  routines  are  based  on  the  follow- 
ing algorithms,   which  express  these  operations  as  rational  algebraic 
combinations  of  the  coefficients. 

Let    P,   Q,     and    R    be  power  series  in  two  variables,    e.  g. 


etc.     (This  notation  was  chosen  so  that  the  formulas  would  be  compati- 
ble with  Fortran.  ) 

If    R    =    P  •  Q,     then  the  coefficient 


■^i 


i  J 

,J     "      Ar^    ^--  ^i,jS-i+l,J"j+l 


If    R  =  P/Q,     then 

I-l       J  J-1 


R 


I' J         S.l 


V.  1-1       j=^l  J-i  > 


where  the  first  sum  is  dropped  if    1  =  1,     and  the  second  if    J  =  1.       If 
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then 


Q    =     -5—  or  R    ^     -r— 

9x  9y 


S.J  =  i-Pm.j    -       ^j  =  ^■''i.j+i 


If 


Q  /       Pdx  or  R 


/ 


then 


\ 


0 


p  •  P 

I"1,J  „  I,J°1 

I  -  1       '  I,  J  J  -  1 


The  series  for    sin  {g(r])-?)    and  cos  {g(n)-?)    were  calculated 
by  considering  the  functional  relations  for    sin  (AtB),     cos  (A+B).      These 
allowed  separation  of  these  series  into  sums  and  products  of  series  of  a 
single  variable.     Other  routines  were  written  to  calculate  the  sine,    cosine, 
exponential,   logarithm  and  power  of  a  series  of  a  single  variable.     These 
routines  were  based  on  recurrence  relations  by  differentiating  products  of 
series.     Leibnitz'  rule  gave  the  formula  to  express  these  derivations. 
Assume-that    P,   Q,   and  R    are  series  in  one  variable,    e.  g. 


^- •      j_'  ' 


i=l 

etc.     (This  notation  was  chosen  so  that  the  formulas  would  be  compatible 
with  Fortram.  )    Suppose 

Q    =    P^ 

then 

Q'P    -    XQP' 

and  '  , .  ' 

dNq'P)    =    XdNqP')      ,  «     .     i-     ,  D^  "^ 


dx     '  ,  k 

dx 
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i=0       ^  i=0       ^ 


so  that  the  coefficients 
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Suppose 

Q    -    log  I  P  I      , 

then 

PQ'     =    P' 

so  that  the  coefficients 


Q^     =    log(P^)     , 
Q2    =     p^     . 


IS.- 1 

|Pk-.i-^Z(i^-)Qk-i.i^ 


Q.   = i^^ (k  =  2, 3, . . . ). 

K+l  _  P^ 

Suppose 

Q    =    sin  (P)  and  R    -    cos  (P) 

then 

Q'     =    P'R  and  R'     =     -  P' Q 

so  that  the  coefficients 

Q^     =    sin  (P^)     ,  R^     =    cos  (P^)     , 

Q2    ^   R.P^     ,  R2    -    -  Q^P^     , 


k-1 

Qia.i      =    RiP,  ^1     +     f    Z^  ^^-^^^■4.^^^     -4-1       ' 
k+l  1    k+l  k    r— r'  1+1     k-i+1 

1  =  1 

k-1 

1  =  1 

If  a  general  function  of  a  function  had  entered  into  the  equations  (5)  or 
(11),    the  problem  would  have  been  much  more  difficult,    especially  since 
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the  coefficients  for  these  functions  must  be  recalculated  each  time 
another  point  of  expansion  is  used.     These  routines  were  important 
because  of  the  form  of    g{0)    in  (12)  and  the  appearance  of    sin  {g(<^)-C) 
in  the  equations  (5)  and  (11). 

5,     The  Body 

After  all  the  flow  quantities  have  been  computed,   the  body  is 
approximated  by  two  methods.     A  stream  function  does  not  exist  for 
this  flow,   so  the  body  is  found  from  the  velocity  field  which  forms  an 
envelope  around  the  body.     The  normal  component  of  velocity  vanishes 
along  the  body.     By  locating  the  stagnation  point  along  the  major  axis, 
the  field  line  that  passes  through  this  point,   perpendicular  to  the  major 
axis,    describes  the  body.     Such  a  point  exists  by  symmetry.     The  field 
line  is  made  up  of  velocity  vectors  with  components    df ,     dn    which  are 
functions  of    u    and    v    only  (see  equation  (7)).     In  practice,   the  field 
cannot  be  found  at  the  stagnation  point,    and  a  starting  point,    displaced 
a  short  distance  (10       radians)  along  the  perpendicular  to  the  major 
axis,   was  used.     This  was  justified  since  the  stagnation  point  along  the 
major  axis  is  assumed  to  be  a  regular  point,    so  that  the  velocity  field 
has  to  be  smooth  there. 

In  the  second  method,   the  body  was  found  by  considering 
A(S)  =  p/p   .      S,     the  entropy,   is  constant  along  a  streamline.     The  body 
is  a  part  of  the  dividing  streamline.     The  body  was  found  by  calculating 
the  entropy  along  a  radial  until  the  value  was  found  to  equal  the  entropy 
on  the  dividing  streamline.     The  entropy  increases  monotonically  from 
the  shock  to  the  body  since  the  shock  monotonically  decreases  in  strength 
from  the  major  axis  to  the  minor  axis.     This  calculation  indicated  that 
the  body  differed  from  the  ellipse 

2  2 

— ^ —  +  — y —  =  1 

(.  32882)^  (.  14226)^ 
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by  less  than    .  006    radians  (see  Figure  4), 

6.     The  Consistency  Checks 

Several  tests  were  used  to  determine  the  accuracy  of  the  com- 
putation.    These  were  also  useful  in  determining  where  to  choose  the 
next  point  for  continuation. 

Bernoulli's  equation  was  not  used  in  the  calculation  of  the 
flow,    so  it  was  available  as  a  consistency  check.     Since  the  Bernoulli 
constant,     B,     is  constant  along  a  streamline,    even  through  a  shock, 
and  since  the  flow  outside  the  shock  is  constant,   then    B    is  constant 

everywhere  in  the  flow.     In  the  initial  conditions    B    was  normalized  to 

2         2  2 

equal  one.     So  by  computing    [27/ (7  -  1  )](p/ p)  +  u     +v     +w      at  differ- 
ent points  and  comparing  this  with    B  -■  1,     a  check  was  provided  on  the 

calculation.     Except  in  the  neighborhood  of  the  body  near  the  minor  axis, 

-7 
the  calculated  value  differed  from    B  =  1    by  at  most    10 

Another  test  consisted  of  comparing  evaluation  of  the  power 
series  expanded  about  different  points  at  common  points  where  both 
series  converged.     As  above,    in  the  neighborhood  of  the  body,    near  the 
minor  axis,    agreement  was  poor. 

The  routine  that  calculated  the  body  from  the  velocity  field  was 
also  used  to  follow  other  particle  paths.     The  calculation  of  the  entropy 
served  as  a  check  on  this  program  and  the  results  were  good  (5 -place 
accuracy)  except  in  the  region  mentioned  above.     In  that  region  there 
are  large  gradients  of  entropy  across  the  streamlines  making  the  test 
unsuitable. 

Significance  arithmetic  was  used  throughout  the  problem.     The 
significance  index  was  a  good  indicator  for  choosing  the  next  point  for 
continuation.     When  evaluating  a  series  the  index  was  examined.     The 
index,    rather  than  gradually  decreasing  as  the  point  of  evaluation  moved 
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closer  to  the  boundary  of  convergence,    would  fall  rapidly  near  the  bound- 
ary;   it  would  stay  constant  inside  the  region  of  convergence. 

Significance  arithmetic  on  New  York  University's  IBM  7090  is 
described  thoroughly  in  [8]  and  [9].     The  ordinary  arithmetic  of  decimal- 
ly represented  numbers  is  carried  out  in  what  is  called  floating  point 
arithmetic;    this  is  without  regard  to  accuracy  or  rounding.     In  the 
machine  a  number  if  represented  by    2   -  f    where    0  <  f  <    1    and    f    has 
27  bitSo     Ordinarily  when  two  numbers  are  added  the  last  step  consists 
of  normalizing  the  answer  by  shifting  the  fraction    f    left  until  the  leading 
bit  is  a  one.     The  exponent  is  decreased  by  one  for  each  place  that    f    is 
shifted.     In  significance  arithmetic  this  step  is  suppressed.     Addition 
and  subtraction  are  followed  by  rounding.     In  multiplication  and  division 
an  effort  is  made  to  give  the  answer  the  same  number  of  zeroes  as  the 
least  significant  operand.     This  is  not  always  possible  due  to  machine 
timing  requirements,   and  so  the  fraction  may  carry  an  extra  leading 
zero.     A  statistical  decision  is  used  in  the  part  of  the  operation  that  ad- 
justs the  number  of  leading  zeroes,    accounting  for  the  discrepancy. 
Rounding  follows  multiplication;,   but  not  after  division,   also  due  to  tim- 
ing.    However  a  trap  follows  division  which  allows  the  programmer  to 
use  a  subroutine  to  perform  the  rounding.     This  was  not  taken  advantage 
of  in  this  program.     The  significance  index  is  27  minus  the  number  of 
leading  zeroes  in  the  fraction,     f.      A  number  with  no  significance  would 
have  an  index  of  zero;    a  nonzero  number  with  complete  significance 
would  have  an  index  of  27. 

7.     Description  of  the  Calculation 

The  problem  was  broken  into  several  subproblems,    each  of 
which  was  solved  by  a  separate  calculation.     The  reason  for  this  was 
that  several  of  the  calculations  took  up  most  of  the  storage  space  avail-   • 
able  in  the  memory;    to  combine  them  would  have  meant  a  loss  in  speed 
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and  flexibility.     Also  many  of  the  tests  were  qualitative,    necessitating 
frequent  interruption  of  the  program.     One  of  these  tests  determined 
the  appropriate  point  for  a  continuation. 

To  start  the  problem,    a  point    (0,  rj    )    and  the  ratio  of  the 
minor  to  major  axis  is  read  into  the  machine  and  the  series  for  the 
initial  curve    g    is  computed.     This  is  an  expansion  in  powers  of    n  -  1^- 
At  this  time    sin  (g(n)-?)    and    cos  (g(r))-?)    are  also  calculated.     These 
series  were  always  calculated  to   24-th  order.     The  coefficients  were 
stored  on  tape  for  the  other  programs.     The  expansions  for  sine  and 
cosine  were  compared  with  the  corresponding  machine  functions  for 
several  different  values.     These  tests  always  showed  agreement  to  at 
least  seven  places. 

The  next  program  calculated  the  expansions  for    p,     p,     u,     v, 
w.      The  order  of  the  series  was  read  into  the  machine  as  well  as 
whether  or  not  the  particular  calculation  was  a  continuation.     The  order 
was  usually  22  at  the  shock,    and  one  less  for  each  continuation.     The 
order  occasionally  had  to  be  reduced  for  expansions  near  the  body  when 
the  machine  had  an  overflow  from  calculating  a  high  order  term.     In  the 
case  of  a  continuation,   the  new  point  of  expansion  had  to  be  given.     Con- 
tinuations were  always  along  radial  lines    (n  =  constant)    through  the 
polar  origin.     It  was  also  necessary  to  recalculate  the  series  for  the 
sine  and  cosine  since  the  equations  were  written  in  terms  of  the  new 

variable    ?'=?-?    .       This  calculation  was  much  simpler  than  the 

o 

original  calculation  for  the  sine  and  cosine,    since  the  functional  rela- 
tions for  these  could  be  used.     The  initial  data  were  written  into  the 
program.     At  this  juncture  of  the  program,   the  operation  for  calculating 
the  series  for  the  flow  variables  is  the  same  as  the  original  expansions 
or  the  continuations.     From  the  Rankine-Hugoniot  conditions  or  by 
evaluating  the  series  to  be  continued,    the  coefficients    p^     .,     p.     ., 
u^     .,     V      .,     w      .    (j  =  1,  ,  .  .  ,  k)    (k  -  1   =  the  order)    are  available. 
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Assume  that  the  coefficients,     p.    .,     p.    .,     u.    .,     v.    .,     w.    .    are  known 

for  all    i  +  i   <   k,     i  <   I.       From  these,    the  coefficients    (dp/drj),    ., 
-  -  '  '         1,  J 

{dp/dr]).    .,     etc.,    of  the  series    dp  I  dr],     etc.,    are  coraputed  for 

^>  J 
i  +  j<k-l,     i<I.       The  coefficients  of  every  series  on  the  right  of 

(5)  are  known  for    i  +  j<k--l,     i<L       Hence    {3p/3C).    .    can  be  com- 

puted  for    i  +  j<k-l,     i<-'-*      This,   then,    gives    p.    .    for    i  +  j   <  k, 

■'-•'J 

i  =  I  +  1     since    idvyl  d^).    .  =  ip.  ,,     ..       Similarly  the  coefficients  for    p, 

1,  J  1+1,  J 

u,     V,     w    are  computed  for    i  +  j<k,     i  =  l  +  l    and  the  process  is  re- 
peated till  all  the  coefficients  are  found. 

After  the  calculation,   these  series  are  stored  on  tape.     This 
program  used  up  most  of  the  memory,   and  could  have  been  written  to 
use  much  less  storage.     However,    as  it  was  coded,    the  program  could 
have  been  converted  to  double  precision  and  use  the  same  amount  of 
storage;    this  was  not  tried. 

Another  program  evaluated  the  flow  quantities,   the  Bernoulli 
constant,    and  the  entropy  along  a  radial  line.     This  was  used,    not  only 
to  compute  tables,    but  to  find  where  continuations  were  necessary,    and 
also  for  the  test  which  compared  power  series  in  overlapping  regions 
of  convergence. 

The  routine  which  calculated  the  streamlines  from  the  velocity 
field  was  based  on  the  following  second  order  scheme. 


Suppose  we  are  at  the  point    (C   ,  n   ).       Calculate  the  point 

o     o 


(?j,  r]^)    by 


=    €_   + 


o  0 


1—2 2         ^ 

/d?     +  dr] 
'       o  o 
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'ii    =   n„   + 
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The  point    (?,  n),     the  end  point,    is  found  by 

dC 
?     =    ?_    + 


°       r^ — 2 


d?.    =    dC 


1 


=    ^o 

1                                                -A 

1 

/      2            2 

ydCj  +dn^ 

d.. 

=   dn 

?=€             (i  =  0,  1); 

n-n. 

d?    and    dr;    are  defined  in  equation  (7);     6    is  a  constant.     This  calcula- 
tion was  done  as  a  separate  program.     The  coordinates  of  an  initial 
point    {^   ,  ri   )    are  read  into  the  machine,    along  with  the  precalculated 
series  coefficients.     At  each  successive  point  the  flow  variables  and  the 
entropy  were  calculated  by  summing  the  appropriate  series.     The  pro- 
gram would  calculate  the  streamline  from    (C    ,  r}   )    to    (C^,  ri    )    using 
expansions  about  a  given    ri  ,     ri     <   ri     <   r\^.      This  program  would 
determine  the  spacing  of  the  expansion  points    (0,  r}   ).       The    r]      were 
chosen  sufficiently  close  together  so  that  the  field  would  be  continuous 
and  smooth  at    (C^,r)   )    and    {^^,r]), 

8.     Results 

Figure  3  illustrates  the  difference  between  the  flows  with 
b/a  =  1,     and    b/a  =  .  7.       Curves  2  and  4  represent  the  shocks,    1  and  3 
the  respective  bodies.     The  particle  paths  for  the  circular  case  are 
shown  in  the  first  quadrant.     Each  streamline  is  a  dividing  streamline. 
In  contrast  with  this  are  the  streamlines  for  the  elliptical  case  shown  in 
the  third  quadrant.     The  streamlines  approach  the  body  fast  and  all  join 
at    0,     which  is  a  singular  point.     In  the  neighborhood  of  the  body  between 
T]  =  1     and    r)  =  7r/  2    the  streamlines  are  bunched  causing  large  gradients 
in  the  entropy;    the  gradients  increased  toward    r)  =  tt/ 2.       This  result 
is  illustrated  in  Figure  5,     The  entropy  is  plotted  against  the  percentage 
of  distance  from  the  shock  to  the  body.     Each  curve  was  found  by 
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computing  the  entropy  along  the  radial  lines    n  =  constant.       The  values 
ri    are  indicated  along  with  the  values    E       of  the  entropy  at  the  shock. 
At    -q  =  l.O    the  effects  of  the  changes  in  entropy  are  apparent.     Addi- 
tional continuations  could  not  improve  the  calculation  near  the  body. 
Between    rj  =  1     and  the  singularity  at    n  =  tt/ 2    the  body  acted  as  if  it 
were  a  line  of  singularities,    and  calculations  were  not  accurate  very 
near  it.     The  expansion  on  the  line    tt/  2    nearest  the  singular  point,    if 
considered  alone,   would  indicate  that  the  singular  point  was  a  stagnation 
point  and  regular.     However  the  velocity  component    v    was  of  opposite 
sign  from    v    calculated  by  neighboring  expansions.     Away  from  the 
neighborhood  of  the  body,    all  the  expansions  were  consistent,   with  high 
accuracy  and  agreement.     As    tt/  2    was  approached  along  the  body,   both 
of  these  properties  were  lost. 

Figure  4  shows  the  body  as  it  was  found  by  calculating  the  en- 
tropy.    It  is  compared  with  the  ellipse: 

2  2 

? +    y .-  1 
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Table  1  illustrates  the  comparison  with  Kopal  [7,    p.    166].     The 
first  column  is  the  argument    0,     given  in  degrees.     The  second  and 
third  columns  show  Kopal' s  calculation  of  the  velocity  components 

(v  =    0).       The  last  four  columns  are  calculations  from  this  problem. 

2  2  7 

B    is  the  Bernoulli  constant    [27/(7  -  l)](p/p  )  +  u     +  w  ,     E  =  p/ p    , 

Table  2  is  a  portion  of  the  tests  used  for  checking  the  power 
series  agreement. 

For  the  circular  and  the  elliptical  case,    a  line  of  singularities 
was  discovered  in  the  region  outside  the  shock  found  by  continuing  the 
functions  expanded  on  the  shock.     This  was  similar  to  the  one  discussed 
by  Van  Dyke  in  [1].       An  atteinpt  at  solving  the  circular  problem  by 
using  a  bilinear  mapping  to  move  the  line  of  singularities  away  was 
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unsuccessful.     A  different  mapping  might  be  more  successful,    but  it 
also  might  be  difficult  to  apply  to  the  elliptical  case;    it  could  introduce 
complicated  functions  into  the  flow  equations. 

Solution  by  power  series  to  the  problem  of  conical  flow  involv- 
ing an  angle  of  attack  should  give  some  interesting  results.     The  singu- 
larity,   found  on  the  body  at    r]  =  it/  2    in  this  problem,    moves  inside  the 
flow  for  the  problem  with  large  yaw.     There  are  also  sonic  lines  in  the 
flow. 
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e 

u  Kopal 

w  Kopal 

u 

15.013 

-. 14140 

.89372 

-.  141398117 

15.0 

-.14102 

. 89376 

-. 141013078 

14.5 

-. 12665 

.89492 

-.126642108 

14.0 

-.  11280 

.89597 

-.112799935 

13.5 

-.09920 

.89690 

-.099197567 

13.0 

-.08567 

.89770 

-.085667565 

12.5  ^ 

-.07210 

.89839 

-.072090499 

12.0 

-.05837 

.89896 

-.0583673328 

11.5 

-.04440 

. 89941 

-.0444039643 

11.0 

-.03011 

.89973 

-.0301010698 

10.5 

-.01535 

.89993 

-.0153441489 

10.  1 

-.00313 

.90000 

-.00311757375 

10.0 

-.00000 

.90000 

.0000090209961 

TABLE  1 

e 

C 

B 

E 

. 37239 

0 

. 99999994 

.0234223269 

0 

.99999994 

.0234223641 

. 36239 

.  0 

.99999994 

.0234347306 

0 

.99999994 

.0234347671 

. 35239 

. 99999994 
.999999985 

.0234513789 
.0234514169 

. 34239 

.99999997 
.999999985 

.0234744802 
.0234745167 

. 33239 

.99999994 
.99999994 

.0235076219 
.0235076576 

. 32439 

1.000000016 
.99999994 

.0235442042 
.0235442936 

TABLE  2A 
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w 


B 


E 


893723622 
993755659 
894922987 
895967498 
896892399 
897699028 
898387432 
898956805 
899405435 
899730816 
899929494 
899994195 
,899996921 


u 


.99999994 
.99999994 
. 99999994 
.99999994 
.999999985 
.99999994 
. 99999994 
.999999985 
. 999999970 
. 999999985 
.99999994 
. 999999985 
.999999911 


0210976347 
0210976191 
0210976191 
0210976250 
0210976243 
0210976250 
0210976243 
0210976236 
0210976236 
0210976236 
0210976236 
,0210976191 
,0210976109 


(Table  1) 


w 


131663367 
131663367 

108251534 
108251497 

085240304 
085240133 

0619579077 
0619578660 

0379564941 
0379568398 

0180792689 
0180794001 


0214699909 
0214702591 

0248421535 
0248424381 

0295664996 
0295668423 

0363407657 
0363413617 

0463652641 
0463660762 

0582596824 
0582522273 


861886069 
861885622 

863083020 
863082573 

864042416 
864041969 

864758521 
864758074 

865212321 
865211874 

865355954 
865355656 


(Table  2A) 
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B 


E 


33626 

0 

.   .999999985 

.0220010728 

0 

.99999994 

.0220010839 

31626 

0 

. 999999985 

.0220816396 

0 

-   .999999985 

-  .0220816530 

29626 

1 

.999999925 

.0222089961 

1 

.99999994 

.0222090080 

27626 

1 

. 99999994 

.0224213980 

1 

.99999994 

. 0224213697 

25626 

3 

.99999994 

.0228058599 

1 

.99999994 

.0228055552 

24026 

3 

.99999726 

.0233712912 

2 

.999999806 

.0233714990 

TABLE  2B 

B 


E 


26529 

0 

.99999994 

.0211160667 

0 

. 99999994 

.0211160667 

24129 

1 

1.00000001 

.0211177625 

1 

. 999999985 

. 0211177647 

21729 

1 

. 99999994 

.0211208947 

1 

.999999970 

■  .0211208962 

19329 

2 

. 99999994 

.0211274773 

2 

1.00000001 

•  ,0211274780 

16929 

2 

. 99999994 

.0211465694 

2 

. 99999994 

.0211466588 

15129 

3 

1.00000001 

.0212243430 

3 

1.00000009 

.0212236516 

TABLE  2C 
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w 


159879729 
159879953 

142547742 
142547995 

127483673 
127479143 

113500647 
113512471 

099161476 
099299692 

086753845 
087303072 


. 0523717880 
.052372162 

.0600173146 
.0600176543 

.067727171 
.067721389 

.074554712 
.074529886 

.079544358 
.079506494 

.082079887 
.082002178 


873486832 
873486608 

876485497 
876485273 

879086256 
879086003 

881223708 
881223485 

882644162 
882645726 

882743970 
882752746 


(Table  2B) 


u 


w 


,143219933 
,143219933 

,115017064 
, 115017317 

,088149212 
,088148281 

,0617579401 
,0617364049 

,0354221612 
, 0350821376 

,0169630229 
,0136737871 


.0112342335 
.0112342350 

.0134353504 
.0134352826 

. 0152707033 
.0152696297 

.0168264821 
.0168320313 

.0178162627 
.0179187737 

.0190556735 
.0167032823 


892936304 
892936304 

896030501 
896030486 

898464277 
898464262 

900253162 
900253013 

901378810 
901375771 

901629582 
901602358 


(Table  2C) 
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B 


E 


26211 

0 
0 

.99999994 
1, 00000001 

.0210980795 
.0210980795 

24011 

1 

1 

,   .99999994 
.99999994 

.0210981175 
.0210981183 

21811 

1 

1 

.999999970 
.999999985 

.0210981823 
.0210981838 

19611 

2 
2 

. 99999994 
.99999994 

.0210983060 
.0210983090 

17411 

2 
2 

. 99999994 
.99999994 

.0210986048 
.0210986122 

15211 

3 

3 

.999999985 
. 999999955 

.0210999608 
.  .0210999362 

TABLE  2D 

Explanation  of  Tables  2 

These  tables  show  a  portion  of  the  tests  on  the  agreement  of 
the  power  series  at  common  points  in  regions  where  both  series  con- 
verge.    The  evaluations  are  along  the  radial  lines    r]  -  constant.       The 
first  column  gives  the  coordinate  along  the  radial  line  in  radians,   the 
next  indicates  which  continuation  of  the  series,       0    indicates  a  series 
which  is  expanded  about  a  point  on  the  shock.     Column  B  is  the  evalua- 


2tP 


tion  of      ,        , , 
(7  -  l)p 

columns  are  the  velocity  components. 


7 


u     +  V     +  w  ,     column  E  is    p/p       and  the  last  three 


Table  2A  compares  expansions  about    rj  =  0,     rj  -  .  2    along  the 
radial  line    n  =  •  1-       The  top  entries  are  from    r)  =  0. 

Table  2B  compares  expansions  about    rj  =  ,  4,     rj  -  ,  6    along 
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u 


w 


■.  141445637 
-.  141445629 

-.114868477 
-.  114868671 

-.089704834 
-.089706264 

-.065109186 
-.065123819 

-.  0405864626 
-.0407251745 

-.0146235406 
-.0160383821 


00179509670 
00179509513 

00213072166 
00213076428 

00241031542 
00241206735 

00264719680 
00266127951 

00279437006 
00289380185 

00316989109 
00261764005 


893703550 
893703550 

896519005 
896519020 

898767680 
898767710 

900469899 
900470048 

901632369 
901633725 

902243361 
902256638 


(Table  2D) 


the  radial  line    r]  =  .  5.      The  top  entries  are  from    r]  =  .  4. 

Table  2C  compares  expansions  about    r7  =  1.3,     r)  =  1.4    along 
the  radial  line    rj  =  1.  35.       The  top  entries  are  from    rj  =  1.3. 

Table  2D  compares  expansions  about    r;  =  1.  5,     r]  -  t:  j  2    along 
the  radial  line    x]  -■  1.5353981629.       The  top  entries  are  from    n  =  1-5. 
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